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ABSTRACT 

An ionization front expanding into a neutral medium can be slowed-down significantly 
by recombinations. In cosmological numerical simulations the recombination rate is 
often computed using a 'clumping factor' that takes into account that not all scales in 
the simulated density field are resolved. Here we demonstrate that using a single value 
of the clumping factor significantly overestimates the recombination rate, and how a 
local estimate of the clumping factor is both easy to compute, and gives significantly 
better numerical convergence. We argue that this lower value of the recombination 
rate is more relevant during the reionization process and hence that the importance 
of recombinations during reionization has been overestimated. 
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1 INTRODUCTION 

Tlie intergalactic medium is liighly ionized by the UV- 
backgrou nd produced by QSOs a n d galaxies, at least since 
z ~ 6 (iGunn fc Peterson! 1 19651 : iHaardt fc Madaul 1 19961 : 
iRauc 3 1 19981 1. Reionization - the transition from neu- 
tral to ionized - occurred around Zreion ~ 10, according 
to the Thomson optical dep th inferred from the cosmic- 
microwave background (see iKomatsu et al.l 120101 . for 7- 
year WMAP result). Reionization starts when the first 
sources of ionizing photons form small, isolated HIT re- 
gions around them. As more and increasingly luminous 
sources form, ionized regions become larger and more nu- 
merous, until they eventually percolate sp ace, signalling the 
end of the epoch of reionization (EoR, lAro ns fc Wingerg 
19721: iHaardt & Madau 1996: Giroux fc Shapird Il996l : 



Gnedin fc Ostrikei] [l997l: IH aiman fc Loebl Il997l. see recent 



reviews by e.g. iBarkana fc LoebT 200ll : ICiardi fc Ferraral 
l20Q5l : lLoeb|[200g ). 

The nature of the sources of ionizing photons is cur- 
rently unknown, with 'first stars', early galaxies, and black 
hole accretion, probab ly all contributing to some extent 
{e.g. iMadau et alll 19991 ) . We recently demonstrated that the 
early population of galaxies predicted by Durham's GAL- 
FORM model produce e nough ionizing photons to complete 
reionization by z ~ 10 jRaicevic et alj|2010al ). The same 
model also matches very well the luminosity function of red- 
shift 2 = 7 — 10 gal axies recently d i scovered by the Hub- 
ble Space Telescope (jBouwens et al.l l2008l : iBouwens et al.l 
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l2009h . In this model, the bulk of photons are produced in 
low-mass (Af, ~ 10® Mq), gas-rich, faint galaxies (rest- 
frame UV magnitude Misoo.ab ~ —16) during a star bursts 
(M* ~ 0.04/i"^Mo yr"^) induced by a merger. 

The fraction TZ of photons emitted per HI necessary to 
complete reionization is 7?, = (1 + A^rec)//eac, where /esc is 
the fraction of ionizing photons that can escape their host 
galaxies and A'rcc is the average number of recombinations 
per HI. We distinguish between recombinations that occur 
in (i) mini-haloes, (ii) Lyman-limit systems (LLS), and (Hi) 
in the general intergalactic medium (IGM). Mini-haloes are 
small high-density clouds that are too cold to form stars via 
atomic line cooling. The presence of mini-halos can have a 
significant effect on the propagation and size of HII regions 
durin g reionization (|Furlanetto fc Ohll2005l : iMcQuinn et al.l 
[2003). However, when overrun by an ionization front their 
;as will be photo-he a ted and they will e ventually evaporat e 
Shapho et all |2004| : llliev et al.ir2005bl : ICiardi et al]|2006l ). 
which decreases their importance in the later stages and af- 
ter reionization. On the other hand, Lyman-limit systems 
have sufficiently high column-densities, A'hi 10^^ cm^^, 
to self-shield. Ionizing photons impinging on such opti- 
cally thick systems get converte d to (non-ionizing) Lyman- 
g radiation at hi gh-efficiency (|Hogan fc WevmannI Il987l : 
I Ranch et al.l |2008| ). These larger systems determine the 
mean fre e path of io nizing photons in the post-reionization 
era (e.g. iMiralda^sc udc 200^. Below we will concentrate 
on the third source of recombinations, those occurring in the 
IGM. 

The recombination rate per unit volume, rtrec, depends 
on the particle density squared, rircc ~ ax^ , where x is 
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the mean ionized fraction, and hence varies thorough out 
the inhomogeneous IGM. Early semi-analytical models of 
reionization used a mean recombination rate, with the IGM 
inhomogeneity expressed in terms of a 'c lumping factor' C, 



hrec ) = cex^ jv?) = a x^C( n)' ^ (e.g. Giroux fc Shapira 



19961: pRgmark et all Il997l: ICiardi fc Ferraral Il997l: 



Hairnan fc Loeblll997l : iMadau et al.lll999l: IValageas fc silM 
19991 ). Numerical simulations of iGnedin fc Ostrikeil (|l997l ) 
yielded high estimates of C ~ 10 (40) at redshifts z = 8 
(5), implying that recombinations are generally quite 
important. Similar values have been used in the estimate of 
the comoving star formation rate density needed to keep the 
post-reionization Universe ionized by balancing ionizations 
with recombinations. The inferred value, 
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{e.g. iMadau et all 1 19991 ) is significantly higher at z ^ 
7 than some curre ntly observationally inferred rates {e.g. 
iBunker et al]|2010l ). 

Current numerical simulations model reionization by di- 
rectly following the propagation of ionization fronts through 
the inhomogeneous IGM, thus in principle eliminating 
the need for a clumping factor {e.g . 
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Ciardi et al 



Trac fc Cen 



2007 



.20031 : llUev et all l2006l : iMcQuinn et al 
I2007D . However the resolution of the radiative 
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transfer (RT) calculation is in general much coarser than 
that of the density fi eld on which the sources are identified 
(e.g. llliev et"aLll2010l ). and a single clumping factor, usually 
derived from higher resolution small box runs, is used to 
take account of the density structure below the re solution 
of the RT mesh (e.q. lCiardi et al.ll2003l : llliev et al.ll2007 ). 

In this letter we will show that a single value of the 
clumping factor is, in fact, not appropriate for estimating 
the recombination rate in numerical RT reionization simu- 
lations, and leads to a significant over estimate of the im- 
portance of recombinations. 



2 DEFINITION OF THE CLUMPING FACTOR 

The recombination rate of Hydrogen-only gas in a volume 
V is 



iVr, 



nldV 



aC{nufV, 



(2) 



where a is the recombination coefficient, nn the hydrogen 
number density, C = {n^) / {nn)"^ the clumping factor, and 
we assumed the gas to be fully ionized (nn = mhii); the 
angular brackets denote a volume average. 

Most current numerical models of reionization follow 
the formation of dark matter structures in a cosmologi- 
cal setting, compute emissivities of galaxies associated with 
dark matter halos, then follow how these galaxies ion- 
ize their surroundi n gs with a radiat i ve transfer calculation 
(ICiardi et all l2003l: llliev et all l2006l: iMcQuinn et all l2007l : 



iTrac fc Cenl2007l : lRaicevic et al.ll2010bl '). Our method in this 
letter is designed to reproduce the steps taken to set up RT 
computational meshes in such simulations. 

We use a dark matter simulation performed with 



the Tree-PM code Gadget-2 (|SpringeJ 120051 '). The simu- 
lation uses 1024=* equal -mass particles in a periodic cos- 
mological volume of size 20 comoving Mpc, assum- 

ing a flat ACDM cosmology with cosmological parameters 
fib, fiA, h, as, Us] = [0.25, 0.045, 0.75, 0.73, 0.9, 1]. 
Baryons are assumed to trace the dark matter, and hence 
the gas density pg is related to the matter density p as 
Pg = {0,b/^lm) p. The matter density p at the positi on Vj 
of particle i is estimated using t he SPH algorithm (|LucvI 
ll977l : lGingold fc Monag"hanlll977^ ■ 
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where m; is the particle mass, and hi its 'resolution length' 
chosen so that it holds ~ 40 neighboring particles j that 
contribute to the sum; W is the smoothing kernel. The (Hy- 
drogen) number density is computed as nn — {1 — Y) pg/rUp, 
where Y is the primordial Helium abundance by mass, and 
TTip is the proton mass. 

Assigning a volume Vi ~ rrii/pi to particle i, allows us 
to compute the mean clumping factor as 
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where A'^part is the number of particles in volume V 
with number density lower than a given threshold den- 
sity nthr- We use rithr = AthA^a) with Athr = 100, 
to exclu de collapsed halos from the IGM density field 



is 



Miralda-Escude et al.l [20001 : iMiralda-Escudd l2003l : 



IPawlik et al.ll2009^ . Baryons in halos do not trace the dark 



matter but collapse to form galaxies. Due to their high den- 
sities, galaxies should not be treated as general IGM, but 
rather as LLS and their effect on the propagation of ioniz- 
ing radiation described in terms of a mean free path, as in 
e.g. iMadau et all (liggTl ). We chose the overdensity thresh- 
old of Athr ~ 100 appropriate for the density at the virial 
radius of a halo, an d also to allow for a direct comparison 
to other works (e.g. IPawlik et al]l2009l ). 

Finally, the RT calculations for large-scale reionization 
models are usually performed on a uniform cubic mesh, with 
the density at each mesh point obtained from the A'-body 
parti cles using, for example , near est grid point interpolation 
(e.g. iHocknev fc EastwoodI Il989l ). We employ several such 
grids in the following discussion. The details of the particular 
simulation we use, and the way we compute densities, are 
not important as far as our conclusions on recombinations 
are concerned. 



3 THE LOCAL CLUMPING 

We will make a distinction between two types of clumping 
factors in the following discussion, both based on Eq. 
The global clumping factor, Cgiobai , is computed by summing 
over all particles, a s is d one in e.g. llliev et al. I (|2005al.l2007f ) 
and iPawlik et all (|2009l ). However we can also divide the 
computational volume in (equal volume, non-overlapping) 
sub- volumes, i.e. a uniform mesl{3, and evaluate C in each 



^ The choice of the volume subdivision is motivated by simplicity. 
Our conclusions are not dependent on the shape of the mesh 
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Figure 1. Evolution of the clumping factor in a simulation with 
1024'' particles in a 20 Mpc box, neglecting particles with 
overdensities higher than a threshold value of Ajj^^ = 100. The 
black solid line shows the global clumping factor, which follows 
reasonably well the result obtained by Iliev et al. (2007; black dot- 
ted line. We then divide the volume in 64'^ equal sub- volumes, and 
compute the local clumping factor in each of them. The dashed 
line shows the mean local clumping factor, with the 50% (99%) 
percentiles indicated by the red (green) shaded region. Clearly 
there is a large scatter in Ciocali and its mean value is signifi- 
cantly lower than that of the global clumping factor. 



sub-volume (mesh cell) by summing only over particles in 
that sub-volum(0. We will call this a local clumping factor, 

^^local • 

We compare Cgiobai and Ciocai computed on 64'^ equal 
sub- volumes of our 1024'^ particles, 20 h^^Mpc aside sim- 
ulation box in Fig. [T] Our value fo r Cgiobai is i n reas onable 
agreement with that obtained by llliev et al.l (|2007l ). even 
though the relations are derived from significantly differ- 
ent N-body runs. Interestingly however, at any z, there is 
a large range of values of Ciocai, up to more than a factor 
of 10 at 2 = 5. In addition also the mean value of Ciocai in 
all sub-volumes of the simulation box is significantly lower 
than that of Cgiobai. Clearly, a single value of C is not able 
to characterize the recombination rate in every region of a 
simulation. 

The mean Ciocai in sub-volumes with a given overden- 
sity A is shown in Fig. (2] This type of C(A) relation was not 
previously discussed in the literature. The highest values 
of Ciocai are found in sub-volumes with intermediate over- 
density, because clumping is a measure of inhomogeneity in 
the density field. Therefore high clumping is found in sub- 
volumes that contain high gradients in the density field, for 
example, small high-density halos in a low density region, or 



^ Note that in general the mean density in each sub-volume will 
be different as well. 
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Figure 2. Mean local clumping factor ((Ciocai)) as a function of 
the overdensity (A) of the sub- volume, at redshift 2 = 5, for vari- 
ous subdivisions of the computational volume; the corresponding 
cell sizes are 625, 312.5, 156.25 and 78.125 /i"^ kpc. Coarser 
sub-cells yield higher values of (Ciocai) at a given overdensity, 
and (Ciocai) peaks at intermediate values of the overdensity. The 
value of Cgiobai is shown as a black dotted line. 



the transition between a filament and a void. Improving the 
sampling of the density field by using smaller sub- volumes, 
decreases the peak amplitude in Ciocai, and moves the maxi- 
mum toward higher overdensity regions, but does not change 
the general shape of the distribution. The general shape is 
also not strongly affected by the choice of Athr. Note that 
the recombination rate, which is oc Ciocai A^, is still highest 
in the highest density regions. However, the use of Ciocai dis- 
tribution shown in Fig. [2] increases the relative contribution 
of moderately overdense regions to the total recombination 
rate. 

An inaccurate estimate of the clumping factor of course 
also implies that the recombination rate is not accurate, and 
hence that the speed of ionisation fronts are not computed 
correctly. Instead of studying the effects on a full RT reion- 
ization run, we show a simpler example by calculating re- 
combination rates in the 20h~^ Mpc, 1024^ particles sim- 
ulation box at redshift z — 5, including densities up to a 
threshold of Athr ~ 100 and assuming that all gas is fully 
ionised. We can sum the recombination rate per particle, 
Nrcc,i = a {m/nip) pg^i, over all particles i, to get the total 
recombination rate, A^rec, correct ~ X^-^rcc.i. This is the 'cor- 
rect' recombination rate as it takes into account all the avail- 
able density field information from the N-body simulation. 
We compare the value of Arcc, correct to values obtained by 
first interpolating particle densities to a mesh as described in 
Section 2, and computing the sum of iVrec in each mesh cell 
using different expressions for the clumping factor (Fig. [3]). 
Not using a clumping factor (yellow bars) shows the effect 
of the density field smoothing by the mesh as the recombi- 
nation rates are both underestimated (a factor of two even 
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Figure 3. Recombination rate, Nrcc, in a fully-ionized simula- 
tion box of 20 h~^Mpc at redsliift 2 = 5 when simulated using 
10243 particles, assuming gas traces dark matter, and imposing 
an overdensity threshold of Athr = 100. The recombination rate 
is computed using various sizes of the sub-volumes over which 
to compute the clumping factor, and is expressed in units of the 
'correct' value computed directly from the A'^-body particles. Not 
using a clumping factor at all (C = 1; yellow bars) leads to an un- 
derestimate of Nrec, which get less at improved sub-sampling of 
the density field. Using a global value (Cgiobai; red bars) leads to 
a large overestimate of N^ec, which gets worse at improved sam- 
pling. Finally using a locally estimated clumping factor (Ciocai; 
blue bars) gives a much more accurate value of Nrcc, which is also 
nearly independent of the sampling resolution. 



at the highest grid resolution) and strongly mesh resolution 
dependent. On the other hand, when the global clumping 
Cgiobai is used to represent the sub-grid matter distribution 
(red bars), it significantly overestimates the recombination 
rate while not remedying the dependence on mesh resolu- 
tion. This is no surprise as a single value of clumping simply 
linearly increases the no clumping result. Crucially, using the 
locally estimated clumping factor, Ciocai (blue bars), leads to 
a much more accurate and resolution independent descrip- 
tion of recombinations (within ~ 25% of iVrec, correct on all 
grid resolutions). Note that the assumption of a fully ion- 
ized simulation box is a special case, chosen for illustrative 
purposes. A more relevant case for the study of reionization 
is a partially ionized cosmological density field. The most 
important consequence of using Ciocai instead of Cgiobai is a 
much lower average recombination rate, for whichever over- 
density regions are ionized at any time. The convergence of 
recombination rates obtained in Fig. [3] with the use of Ciocai 
also leads to the converge nce of I-front speeds du ring reion- 
ization as we will show in I Raicevic et al.f()2010bl ). where we 
also discuss the ionized fraction as a function of overden- 
sity during reionization. Also, Ciocai is computed assuming a 
full ionization and does not provide a perfectly accurate re- 
combination rate estimate for partially ionized sub- volumes. 
Therefore, the use of a higher resolution RT computational 



mesh is always preferable to employing any kind of clump- 
ing. 

4 CONCLUSIONS 

Clumping factors are often used in simulations of reion- 
ization to represent the unresolved matter inhomogene- 
ity, below the computational mesh resolution, which may 
significantly contr i bute to the recombination rates (e.g . 
ICiardi et aLll2003l : lliiev et al.ll2007l : iMcQuinn et al.ll2007l ). 
Here we have shown that, because the density field during 
reionization is so inhomogeneous, using a single clumping 
factor in general leads to a significant over estimate of the 
recombination rate (factors of several), which may in fact get 
worse as the grid resolution is improved. As a consequence 
the speed of ionization fronts is artificially depressed, reion- 
ization delayed, and the outcome of the radiative transfer 
sim ulations significantly re solution dependent. 

iMcQuinn et~all (|2007l ) improved on this issue by ana- 
lytically deriving the clumpi ng factor as a functi on of over- 
density, C = C(A), and also iKohler et all (|2007l ) took into 
account some density dependence by splitting their simula- 
tion volumes into 8 equal sub-volumes. However, both ig- 
nored the fact that clumping still depends on the volume 
over which it is computed {i.e. the size of the sub- volumes), 
and that dependence is not negligible as shown in Fig.[2l Not 
taking the volume into account will always lead to the RT 
results depending on the computational mesh resolution. 

Fortunately it is not computationally intensive to com- 
pute a local clumping factor on a reasonably fine grid. The 
recombination rate obtained using this clumping factor is 
very close to the 'correct' value inferred directly form the 
particles themselves, and is not very sensitive to the resolu- 
tion of the mesh used (see Fig. [3|. In our dark matter only 
simulations the cell size should be such that most of the 
structures are resolved, and should id eally be close to th e 
Jeans mass of the photo- ionised IGM (jPawlik et al.ll2009l ). 
We therefore argue that any future fits of clumping factors 
aimed for use in RT simulations of reionization must in- 
clude both density and volume dependence. The alternative 
of computing the recombination rate per particle directly 
(jTrac &i Cenll2007l ). is far more computationally expensive. 
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